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Abstract 

The evolution of a system of chemical reactions can be studied, in the eikonal approximation, by 
means of a Hamiltonian dynamical system. The fixed points of this dynamical system represent 
the different states in which the chemical system can be found, and the connections among them 
represent instantons or optimal paths linking these states. We study the relation between the 
phase portrait of the Hamiltonian system representing a set of chemical reactions with constant 
rates and the corresponding system when these rates vary in time. We show that the topology of 
the phase space is robust for small time-dependent perturbations in concrete examples and state 
general results when possible. This robustness allows us to apply some of the conclusions on the 
qualitative behavior of the autonomous system to the time-dependent situation. 

PACS numbers: 05.40.-a, 05.45.-a, 64.60.My, 82.20.-w 
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I. INTRODUCTION 



Understanding the dynamics of chemical kinetics is a question of fundamental character 
in nonequilibrium statistical mechanics, apart from being of broad importance in applica- 
tions to other sciences. Indeed, many models in chemistry jl|, biochemistry j^j], ecology 
and biology j|| use stoichiometric relations as theoretical first principles describing some 
phenomenon. The simplest description of the time evolution of a set of N reacting species 
is probably given by the mean-field equations, an N— dimensional dynamical system repre- 
senting the concentrations, or total number of molecules of the reacting species. One of the 
advantages of this approach is that it allows us to use the powerful machinery of dynamical 
systems theory 5|, raising the possibility of identifying stationary states with fixed points, 
periodic behavior with limit cycles, etc. 

Of course, as with every theory in physics, the mean-field approximation has a range of 
validity. As it ignores fluctuations, its description of the chemical system might be accurate 
for short times; however, long-time dynamics will be affected, dramatically in some cases, by 
rare events. It is possible to study exactly the system evolution as a time-continuous Markov 
process, the probability distribution of which is given by the solution of an adequate master 
equation f], Q] . The full analytical solution of the master equation is, in most situations, 
a formidable problem, and it yields, usually, much more information than that needed in 
applications. The development of approximate theories is thus clearly justified. 

One of the most common approximations is the Fokker-Planck equation, obtained from 

n q 

the master equation by means of a Kramers- Moyal or van Kampen size expansion |6|, |7|. This 
theory assumes that the number of reagents is very large, so we can consider the implicit 
stochasticity of the process as small Gaussian fluctuations around the mean-field behavior. 
Definitely, rare events do not belong to this category. If one wants to deal with fluctua- 
tions that are comparable with, or even greater than, the mean value, any approximation 
scheme must not rely on assumptions such as an asymptotically large value of the number of 
reagents. It is possible to construct such a theory, for instance, by taking advantage of the 
relation among chemical kinetics and quantum mechanics 8|, |9[ . Once the master equation 



is formulated as a quantum problem, it is possible to develop eikonal 10, 111 or WKB ap- 

□ n 

proximations [12j, or matched asymptotic expansions of the spectral formulation [13j, able 
to tackle rare events. In this work we will concentrate on the eikonal approximation intro- 
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duced in 



111 ]. One of the advantages of this approach is that it reduces the problem again 



to a dynamical system of 2N dimensions for N reacting species. The additional N degrees 
of freedom are the conjugate "momenta" corresponding to each of the concentrations and 
represent a measure of the size of the fluctuations. Due to the Hamiltonian symmetry of this 
system, it can be effectively reduced to a (2N — 1)— dynamical system on some Riemannian 
manifold. So, if there is only one chemical species, the case in which we will concentrate 
here, we have a reduced number of dynamical scenarios. In particular, only fixed points will 
be of physical relevance, and they will denote the possible stationary states in which the 
system can be found. Contrary to the mean-field situation, rare events can drive the system 
for one stationary state to another, a fact that is reflected by the existence of connections 
between the fixed points of the dynamical system. These connections are not the unique, 
but the optimal way in which a system evolves from one state to the other [lo|, and, as they 
are reminiscent of instantons in quantum mechanics [ljj], we will denote them as instanton 
connections. Its importance is huge: the web of connections, having the fixed points at 
their intersections, encodes the qualitative behavior of the chemical system H], and its 



topology serves as a princip 



reaction-diffusion models |15] 



e for the classification of nonequilibrium phase transitions in 



All these works have been focused on the dynamics of chemical kinetics happening at 
constant rates. While this assumption is reasonable in many cases, there are situations 
in which we should go beyond it and consider the explicit time variation of the reaction 
rates. Periodically illuminated chemical reactions or seasonal variation in population 
dynamics [l7] serve as examples of nontrivial behavior generated by a temporal forcing. 
Not much attention seems to have been paid to eikonal approximations of time- dependent 
chemical kinetics, and when they are considered, nonautonomous perturbations are usually 
treated within the quasistationary approximation Q| • It is our goal to extend the exis- 
tent approaches and consider arbitrary frequency, albeit small, perturbations. The concrete 
problem under study is how the phase portrait of the eikonal dynamical system is modified 
when time- dependent perturbations enter into play: do the instanton connections survive 
or do they disappear changing the qualitative behavior of the system? Giving a totally gen- 
eral answer to this question is a difficult task, but we will see that these connections seem 
to be very robust and persistent when they are subject to a small periodic forcing. This 
will be shown with the help of particular models, for which rigorous results (for the eikonal 
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FIG. 1: Phase space for the branching and annihilating particle problem. The parameter values 
are a = X = 1. 

approximation) are easily proven, and then we will extend them to the general setting when 
possible. 

II. INSTANTON PERSISTENCE 
A. Branching and annihilation 



For the shake of clarity, we will illustrate the problem of instanton persistence with 
a particular reaction, branching and annihilation of identical particles, but we will state 
general results for arbitrary reaction sets [20f] . Let us start considering a single species of 
identical particles A, which annihilate in pairs and undergo binary branching 



A + A 



A^A + A. 



(1) 



The master equation describing the probability distribution of having n reagents at time t 
reads 

= a(t)[(n - l)P n -!(t) - nP n (t)} + *@-[(n + 2)(n + l)P„ +2 (t) - n(n - l)P„(t)]. (2) 

It can be represented by a partial differential equation (PDE) by introducing the generating 
function 

oo 

G{p,t)=J2p n Pn(t), (3) 
n=0 



which provides us with the time-dependent Hamiltonian 

H(p, q, t) = a(t)(p - l)pq + ^(1 - p 2 )q 2 (4) 
and the imaginary time Schrodinger equation 

d t G=-H(p,-d p ,t)G. (5) 



The eikonal approximation proposes the reduction of the problem to Hamilton equations , 
which in this case become the nonautonomous dynamical system (note, however, that 
we could reduce the problem to one time- dependent reaction rate by means of the time 
reparametrization t = f o~(r)dT) 

OH 

P = —^ = v(t)(l-p)p + \(t)(p 2 -l)q, (6) 
q = ^ = cr(t)(2p-l)q~X(t)pq 2 . (7) 

Suppose for a moment that both a and A are time independent. In this case the system 
exhibits three lines with zero energy: the invariant lines p = 1, q — 0, and 

Furthermore, we know that these three lines connect the fixed points (0,0), (1,0), and 
(l,<j/A) in the (p, q) plane, see Fig.([T]) (or Fig.(3) in [11]). Now we can try to understand 
what happens when we let the reaction rates vary in time. It is easy to see that the lines 
p = 1 and q = remain invariant zero-energy lines of the system; however, the explicit time 
dependence of the Hamiltonian prevents the conservation of the energy H and, in general, 
forces the disappearance of the third zero-energy line. We can generalize this fact for an 
arbitrary reaction Hamiltonian. Given any set of reactions with time-dependent rates, we 
know that the Hamiltonian necessarily fulfills 

H(p=l,q,t)=0 (9) 



due to the conservation of probability III, |l5|]. So this means that for any set of (time 



dependent) reaction rules, the p = 1 line is an invariant, zero-energy line of the dynamical 
system. Since this line describes the mean-field dynamics of the system, we will call it the 

nn 

mean-field line [111 . 1191 ] . Some systems possess an absorbing state when they contain zero 



particles; this happens if all reactions of the type 



nA (10) 



are absent in the dynamics. In this case, the Hamiltonian must obey the condition ll|, [l5j 



H(jp,q = 0,t) = 0. (11) 

So we can claim that any (nonautonomous) system without particle production from the 
vacuum has the invariant, zero-energy line q = 0. And thus, when it is present, we will call 
it the absorbing-state line. These properties can be clearly seen from the general form of the 
Hamiltonian term representing the reaction 

mA nA . (12) 



it is [1 



H mn = ^^(p n -p m )q m - (13) 
ml 



In the autonomous situation, the set of zero-energy lines determines the physics of the prob- 
lem: the fixed points obtained when these lines cross represent the possible states in which 
the system can be found and the connections among them the possible transition paths. 
Together with the mean- field line (the global minimum of the action) and the absorb ing- 
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19|]. We know 



state line one finds other lines with zero energy: the instanton lines 
that both the mean-field line and the absorbing-state line persist if we let the reaction rates 
vary in time, but however, it is not so obvious to see what happens with the instanton lines. 
These are defined in terms of zero energy if the system is not explicitly dependent on time, 
but when this is not the case, the definition loses its meaning since energy is no longer con- 
served. Due to this fact and because the physical role of the instanton lines is to be optimal 
paths between different states, what we would like to know at this point is if both fixed 
points and connections among them survive after the nonautonomous forcing is switched 
on. We can be sure that hyperbolic fixed points persist to a small periodic time-dependent 
forcing, after possible relocation of their position, but do the instanton connections persist? 

To address this question we will use the method developed by Melnikov 21] (see also Q]). 
Consider an autonomous two-dimensional dynamical system 

x = f(x), (14) 

x = (x\, X2), f = (/1, ^2), with two hyperbolic fixed points x a and x?, linked by a heteroclinic 
connection, which is parametrized by the system solution x/j(t — t ) f°r initial time to- 
Assuming that the motion goes from x a to x?,, this connection is formed by a branch of the 




FIG. 2: Sketch of the phase space for the branching and annihilating particle problem, assuming 
that the stable and unstable branches intersect. The heteroclinic connection goes through all the 
intersections of these branches. 

unstable manifold of x a , which totally overlaps with a branch of the stable manifold of x&. 
Let us now consider the perturbed version of this problem, 

x = f(x)+eg(x,t), (15) 

where the perturbation g(x, t) is time periodic with period T, amplitude e small enough and 
sufficiently regular. The dynamics of this nonautonomous system is given by the associated 
Poincare map V e , which maps every initial condition point x(0) with the corresponding 
value of the solution x(T) after one period has elapsed 5j. The hyperbolic fixed points of 
the unperturbed system, x a and Xf,, are hyperbolic fixed points of the Poincare map Vq. 
Since the map V e is a perturbation of Vo, these points have a continuation, x^ and x£, as 
hyperbolic fixed points of V e , and their invariant manifolds vary continuously with respect to 
e. The perturbed system will not, in general, maintain the coincidence between the branches 
of the unstable and stable manifolds of x a and x&, respectively: now these branches might 
intersect, preserving the existence of the heteroclinic connection (see Fig.(j2J)) or might not, 
destroying it, like in Figs.® and (j4]). The distance between the stable and unstable branches 
is given by d(e,t ) = eM(t ) + 0(e 2 ), with 

/oo 
f (x h (* - t )) A g(x ft (t - t ),t)dt, (16) 
-oo 
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FIG. 3: Sketch of the phase space for the branching and annihilating particle problem, assuming 
that the stable and unstable branches do not intersect. 



where 



f Ag 



fx h 
9i 92 



(17) 



denotes the wedge product of vectors f and g. Equation (|T6|) defines the so-called Melnikov 
function, which yields the first-order approximation in e of the distance between the stable 
and unstable manifolds measured along a direction that is perpendicular to the unperturbed 
connection at the point X/j(£q). A change of sign of M(to) means that there exists some 
to such that d(e, to) = 0, implying the existence of a solution x. e h (t) of Eq. (JT5l) denning a 
heteroclinic connection among the two hyperbolic fixed points x. e a and x^ of the Poincare 
map corresponding to Eq. (fl5j) . say, 



lim x£(t) = x* 



lim x£(t) = x£ 



(18) 



Let us now consider the branching and annihilating system, Eqs. ([6]) and (J7J), with con- 
stant reaction rates. The instanton connection (the separatrix linking (l,o~/A) to (0,0)) is 
parametrized by the system solution 

1 2 ^ X ^ (19) 

(20) 
(21) 



x fc (t - to) = (Ph(t - to), q h (t - t )) - x i + eCT(t _ to) , 2 + eCT(t _ to) 
If we assume that the perturbation of the reaction rates is given by 



a(t) 
X(t) 



o + eax(t), 
A + eAi(t), 
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FIG. 4: Sketch of the phase space for the branching and annihilating particle problem, assuming 
that the stable and unstable branches do not intersect. 



and the new rates are —periodic and regular enough, we obtain the system 

p = a (l- p ) p + X(p 2 -l)q + e[a 1 (t)(l-p)p + X 1 (t)(p 2 -l)q] 
q = a(2p - l)q - Xpq 2 + e [a 1 (t)(2p -l)q- Xi(t)pq 2 ] , 



(22) 
(23) 



which is of type (1151) . The associated Poincare map V e still has the origin as hyperbolic fixed 
point. Another hyperbolic fixed point of this map, corresponding to (1, er/A) when e = 0, is 
(1, q e ), where q e is obtained through the solution of Eq. (123!) for p = 1. This is nothing but 
a Bernoulli differential equation, which can be straightforwardly integrated to get 



q(0) exp J* a(r)dr 



1 + q(0) Jo exp [J Q S a(r)dr] X(s)ds 
We just need to apply the condition q(0) = q(T) to this solution to find 



(24) 



exp 



Jo a ( T ) dT 



- 1 



Jo ex P U S(r ( T ) dT } X ( s ) ds 



(25) 



To determine if the unstable manifold of (1, q e ) intersects the stable manifold of (0, 0) we 
substitute the corresponding values of f and g into the Melnikov function (|T6j) . 



M(t ) 



[A<7 X (t) - a\ 1 (t)] [(1 -p)(l-p- p 2 )q 2 ] ( Xft (t - t ))dt = 

POO 

h(t)(j)(t - t )dt = / h{s/a + t )4>{s)ds, 



(26) 



after the change of variables s = a(t — to) and where 



h{t) = Aoi(t) -crAi(t), 

<f>(t) = [(l-p)(l-p-p 2 )q 2 ](x h (t)), 
~ 4a 2 e 2s [l + 2sh(s)] 



rw A 2 (l + e s ) 3 (2 + e s ) 2 ' 

the symbol sh standing for the hyperbolic sine. It is easy to see that </>(s) has zero mean 

e 2s [l + 2sh(s)] 



(27) 
(28) 
(29) 



[s)ds 



4V 



:ds 



2a 2 
"A2 



_(1 + e s ) 2 (2 + e s )_ 



0. 



(30) 



, (1 + e s ) 3 (2 + e s ) 2 

a fact that will prove its usefulness later on. Since h(t) is T periodic and continuous, we can 
expand it in Fourier series, 



E 



2-Kint/T 



h(t) = > r;„f 

n=— oo 

and substitute it into the Melnikov function to obtain 



(31) 



M(t ) 



4V 
" A 2 



E 



2nint Q /T 



oo ^2s 



e 2s [l + 2sh(s)] 
'1 + e s ) 3 (2 + e s ) 



(32) 



which is the Fourier decomposition of the Melnikov function. We can check the validity of 
the Fourier series (|32l) after contrasting that the norm 



e 2s [l + 2sh(s)] 



:i + e*) 3 (2 + 



ds 



5^5-11 



(33) 



is finite. This representation allows us to calculate the mean value 



T 



M(t )dt = a / (f)(s)ds = 



where we have used 



T 



e 2mnt /T dto = $ 



(34) 



(35) 



and 5 n o denotes the Kronecker delta. So we know that the Melnikov function is periodic, 
continuous (from its very definition in Eq. ffl6l) it only depends on f, g, and X/J, and with 
zero mean, implying that it either crosses zero or vanish identically. In the first case, we 
can claim that the point (1, q e ) is connected to the origin, in the second, we cannot, due to 
the possible presence of small terms 0(e 2 ). However, we can see that the second case only 
happens when h(t) is constant. Indeed, integrating by parts and changing variables x = e s 
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we obtain 



oo 9s 



e 2s [l + 2sh(»] 
'l + e s ) 3 (2 + e s ) 



2nins/(aT)^ s 



nxn 



aT ./_«, (l + e s ) 2 (2 + e* 



-e 2 ™ s/((TT) cis 



aT 7 (l + x) 2 (2 + a; 



-cix. (36) 



This last integral can be obtained by means of contour integration in the complex plane. 
Noticing that for a complex variable z 



2mn/(aT) _ 2nin\n(z) / (aT) 

A/ O « 



(37) 



we can use the keyhole contour choosing the logarithm branch cut as the positive real axis, 
and by employing the residue theorem we obtain 



x 2irin/(aT) 

o (1 +x) 2 (2 + x) 



dx 



2m 



{Res[/(^) ) -1] + Res[f(z), -2]}, (38) 



where 



I _ e -4vr2n/(o-T) 

„2mn/(aT) 



We can now compute the residues 



Res[/(,z), — 2] = exp 
Res[/(z),-l] 



2ix 2 n 2irin 
+ 



aT 

^ z 2mn/{oT) 



aT 



dz (2 + z) 



ln(2) 



1 + 



2mn 
~~o~T 



exp 



2it 2 n 
' aT 



(39) 

(40) 
(41) 



to conclude 



2n 2 n/(aT) 



exp 



27r 2 n 



00 e 2s [l + 2sh(g)] 
_ (i + e s)3 (2 + eS )2 



2nin 
1 H — exp 



e 2nins/(aT) ds 

2-ftin 



aT 



ln(2) 



(42) 



1 — exp[— 4ir 2 n/(aT)] r ^ aT J \ aT 
and one can see that it only vanishes if n = or n —>■ ±oo, for a and T positive finite real 
numbers. Hence, the Melnikov function is identically zero if and only if a n = for all n ^ 
in Eq. ( 13~TT) or, what is the same, if and only if h(t) is constant. 
We still have to analyze what is the meaning of the condition 



Acxi(t) — a\i(t) = c, 



(43) 



for some constant c. If c = 0, then the system can be reduced to the unperturbed one by 
means of a reparametrization of the time variable: 





i+ u ; 




Jo 


a 



dr. 



(44) 
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In this case the geometry of the phase space is exactly preserved; only the parametrization 
of the system solution along the separatrices might change. If <j\ and Ai are chosen constant, 
then for any value of c, the phase-space topology is the same, as a small retuning of the 
Hamiltonian parameters cannot modify it. In these two particular cases we know that the 
Melnikov function is identically zero because the connection is still given by the complete 
superposition of a branch of the stable manifold of one of the fixed points and a branch of 
the unstable manifold of the other. In the general case we know that a perturbation fulfilling 
hit) = c is a combination of these two, translation and reparameterization, 

a(t) = (<7 + e<7i)[l + e0(t)], (45) 
A(t) = (A + eAi)[l+e0(t)], (46) 

for some function <p(t), up to terms 0(e 2 ). This means that, in order to know if the distance 
between the invariant manifolds of the fixed points is identically zero, we would have to use 



the extension of the Melnikov method to higher orders 23]. 



B. General setting 

It would be highly desirable to extend two properties of this example to general systems, 
say, the possibility of expressing the Melnikov function as a Fourier series, and to show that 
it has zero mean. We will start with the second claim assuming that the first one is true, 
and we will check its validity afterwards. The most general reaction Hamiltonian can be 
built adding the generic terms (TT31) . 

H=J2M^(P n -P m )q m , m,n = 0,1,2,-.., (47) 
^— ' ml 

and it allows us to express the dynamical system as 

J2^%( P m -p n )Q m -\ m,n = 0,1,2,-.., (48) 
*—f (m — 1)1 

m^pn 

^ = Am^*)/ „-l _ m-1) m m,n = 0,1,2, ■■■ , (49) 
ml 

where all the rates (3 mn {t) are supposed to have the same period T (the case of perturbations 
with different periods will be discussed below). Assuming the form of the perturbation 
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Pmnif) = Pmn + e Pmn(t) allows us to split the Hamiltoiiian into two terms H = Ho + eHi- 
So we can write the integrand of the Melnikov function as the set of Poisson brackets 

f(xh(*-* )) A gfoft -*„),*) = {H!{yi h {t-t Q ),t),H Q {^ h {t-t ))} (50) 

and obtain 

/oo 
{Wi(xfc(t - *o), *), Wo(x ft (i - to))}dt = 
-oo 

J —Hi(x h (t-t ),t)dt- j —Hx{x h {t-t G ),t)dt = 



-oo 

"OO 



/•OO Q 

J —Hx{x h {t-t ),t)dt, (51) 



where we have used the fact that instanton lines link fixed points lying on zero-energy 
invariant lines. Since 



Wl(Xfc(*-*o), *) = J2 Pmn(t)Jmn(t-t ), J m n(f-t ) = -^(p(t-t ) n -p(t-t ) m )q(t-t ) m , 

*—f ml 

(52) 

we can integrate by parts to get 

Y<J t P™{t)3rnn{t-t )dt= J ^(t) - ^(t - t )^. (53) 

This allows us to rewrite, at least formally, the Melnikov function as a Fourier series, 

°° POO J 

M(t ) = £ e2mkt0/T J T,^j-Jrnn(s)e^ ks/T ds, (54) 

after the change of the integration variable s = t — t . Its mean value is thus given by 

1 f 00 d 

- I M{t )dt = I Plnfi^Jrnn&ds = 0. (55) 

In order to check the validity of Fourier series ( 1341 we need to bound the norm 

dt. (56) 

This is not a difficult task because both p(t) and q(t) must be continuously differentiable 
along the connection and with finite limits when t — > ±00, which must coincide with their 
values at the fixed points. Furthermore, the derivatives p and q decrease to zero exponentially 



/■oo 


d^Tmn 


/ —00 


dt 
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at infinity due to the hyperbolic character of the fixed points. This implies that the norm (|56|) 
is necessarily finite. 

It is actually simpler to show that the Melnikov function has zero mean. Equation ([55]) 
can be cast into the form 



M(t ) 



-IT 22P™(t + to)Jmn(t)dt, (57) 

ato J-oo , 



which, together with the periodicity of the functions /3^ n (i), directly provides the desired 
result. However, the Fourier series setting favors the identification of the conditions un- 
der which the Melnikov function vanishes, as Eq. (143]) in the last section, which cannot 
be straightforwardly extracted from this expression. So we believe that the Fourier series 
representation will be more helpful in applications to concrete reaction schemes. 

In conclusion, we have shown that the Melnikov function of a perturbed instanton con- 
nection can always be expanded in Fourier series, provided that it links hyperbolic fixed 
points. Furthermore, it is always a zero mean function, which implies that it either crosses 
zero or vanish identically. So for any reaction Hamiltonian we have the following result: 
given an instanton connection linking two hyperbolic fixed points, any small time-periodic 
perturbation of the rates will preserve the existence of the connection if the corresponding 
Melnikov function is not identically zero. In the other case, we would have to study the 
behavior of the Melnikov distance at higher orders, in order to be able to give a rigorous 



conclusion 



23]. 



To finish, let us note the difficulties in extending this result and proving a general the- 
orem about the persistence of the instanton connections even when the Melnikov function 
is identically zero. In this case, we could recall the expansion of the Melnikov distance in 
terms of higher-order Melnikov functions [3| . One would be tempted to use an equivalent 
argument to that of the present section to try to show that an arbitrary-order Melnikov 
function either has zero mean or vanish identically. This would imply in turn that either 
some function in the expansion is nonzero, and thus the connection is preserved by means 
of a transversal crossing of the stable and unstable manifolds, or the whole series becomes 
identically zero. While this can suggest, at first sight, that the Melnikov distance is also 
identically zero in this case, we cannot rely on a rigorous argument to prove so, since the 
perturbative expansion is not analytic in the small parameter in general situations. Fur- 
thermore, we can always find a family of transformations (for instance, by continuing the 

14 



chain of reparametrizations in Eqs. (145!) and (146]) ) able to nullify an arbitrary order Melnikov 
function. So a full proof would need to handle these special perturbations, and this is com- 
plicated by the lack of analyticity of the series expansion, despite the probable fact that they 
are built by combinations of reparametrizations, constant shifts of the parameter values, or 
other trivial transformations preserving the phase-space topology. Also, proving that the 
higher-order Melnikov functions have zero mean is not a trivial fact. These imply nonlinear 
combinations of the external perturbations and the corresponding mixture of Fourier modes 
that complicates the development of a simple form like Eq. (1541 in these cases. We illustrate 
this situation with the second-order Melnikov function in the Appendix. Anyway, perturba- 
tions requiring a higher-order Melnikov analysis do not constitute the generic case, and we 
believe that the majority of the physical situations could be analyzed within the first-order, 
or at most second-order, formalism. 



III. DISCUSSION 



So far we have concentrated in showing that the instanton connections persist when the 
chemical system is temporally forced by a weak perturbation. We devote this section to 
explain the physical consequences of this fact. As we already argued, the topology of the 
phase space describes the physics of the system, so if its topology is preserved when the 
system is forced, this means that the qualitative behavior of the system is preserved too. 
The disappearance of an instanton connection would mean that this qualitative behavior is 
modified, but how strongly? These connections are optimal paths linking different physical 
states (given by the fixed points of the dynamical system), but the absence of one connection 
communicating two states does not mean that the system cannot go from one to the other. 
It only means that there is not an optimal way of going. We will show the importance of 
this fact using a toy model of biological relevance in plankton modeling. Suppose we have 
the reactions 

A^UA + A, (58) 
occurring at the same constant rate 7. This set of equations has been used to model plankton 



patchiness in some occasions 
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equation for the generating function 



26 L We can straightforwardly write the Schrodinger 
270, 



d t G = 1 { P -ifd p G, 



(59) 
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so the reaction Hamiltonian reads 



H = j(p — l) 2 g. (60) 

The corresponding dynamical system is highly degenerate, with two invariant lines p = 1 and 
q = 0, with all their points being fixed. This means that if we start with an initial condition 
(q = no,p = 1), for some n > 0, then this will be the solution for all times. As is totally 
clear, there are no instanton lines linking the mean-field line with the absorbing-state line. 
Let us now concentrate directly in the exact PDE fl59l) instead of the approximate dynamical 
system. This is a linear, first-order PDE which can be easily solved along characteristics [6J] 



G(p,t) = G 



' l-p 



l + (l-p)jt 

and due to normalization, we have the long-time asymptotics 

1 — p 



lim G(p,t) = lim Go 



1 - 



(61) 



G (l) = 1, (62) 



l + (l-p) 7 t_ 

which corresponds exactly to the probability distribution P n = 5 n Q. This means that in the 
infinite-time limit, regardless of the initial condition, the system will be in the absorbing 
state with zero particles. So, as we have seen, the absence of instanton connections does 
not stop the system from undergoing an extinction transition. One can figure out how this 



happens representing the process fl58|) with the Ito stochastic differential equation 



dp = ^fpdW, (63) 

where W denotes a Wiener process. This equation describes Brownian motion with state- 
dependent diffusion and tells us that the transition to extinction is very different in this 
case: no optimal trajectory is chosen; instead, the state with zero particles is reached after 
performing a random walk from the initial condition. This will not necessarily happen in ev- 
ery situation; actually, state-dependent Brownian motion is one of the simplest possibilities. 
More complicated Hamiltonians with powers of the momentum above the second exhibit 
properties corresponding to non-Gaussian statistics [(J. In general, the random walk will 
be more complex than simple Brownian motion, but all the cases will have in common the 
absence of an optimal way of going from one physical state to another. 

The instanton connections are a useful tool that allows us to calculate the frequency of 
rare fluctuations that drive the system among different states. The mean transition time is 

16 



obtained, at exponential order, computing the action along the instanton connection 11]. 
If our perturbation is pure time reparametrization, as the second perturbation in Eqs. (|45p 
and f)46p . and the function periodic, a straightforward computation reveals that the 
action along the heteroclinic orbit is exactly the same as for the unperturbed system. More 
general perturbations are not tractable analytically, and we would have to rely on a numeri- 
cal treatment of the problem. The heteroclinic connection can be obtained, for instance, by 
means of a shooting method and the action computed as a numerical integral on it. In any 
case, we expect that small, well-behaved, periodic perturbations will yield transition times 
of the same order of magnitude as the unperturbed system. Instanton persistence is impor- 
tant as it allows approximating this sort of nonautonomous systems by their autonomous 
counterparts in the first instance and serves as a justification of the quasistationary approx- 
imation for slow signals. Perturbations that are not periodic might have a stronger effect on 
the system, and in particular, singular enough perturbations could change the phase-space 
topology; however, it is rather difficult to conceive physical situations that give rise to such 
a perturbation. Also, the mathematical framework changes considerably, since the absence 
of a Poincare map does not allow the definition of the instanton as a connection among fixed 
points of this map, like in the periodic case. In any case, even if the instanton connections 
were broken, but the physical states did not drift apart form each other, transitions will be 
possible if we wait long enough, due to the uncorrclatcd fluctuations the system is subject 
to 28j, but however, there will be no optimal paths linking those states. 



IV. CONCLUSIONS AND OUTLOOK 



In this work we have studied the persistence of instanton connections in chemical systems, 
when we promote the reaction rates from constants to functions of time. A set of chemical 
reactions can be mapped, under the eikonal approximation, onto a dynamical system, the 
fixed points of which denote the possible states where the chemical system can be found. 
Connections among fixed points denote optimal transition paths communicating different 
states. The persistence of these instanton lines indicates the robustsness of the physical 
processes that are taking place in the reactor and shows us that the qualitative behavior of 
the system will be similar when subject to small nonautonomous perturbations. We have 
shown our results with one particular model, and we have stated them in a more general 
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setting when it has been possible. The main theoretical tool that we have employed is the 
Melnikov function, which has proven itself very useful in dealing with this type of problems. 

One of the directions in which we would like to extend the theory is the problematic of 
connections linking nonhyperbolic, or one nonhyperbolic and one hyperbolic, fixed points. 
The physical motivation comes from the recent classification of phase transitions in reaction- 
diffusion models, which identifies critical points of the physical system with the nonhyper- 



bolic fixed points at a bifurcation threshold of the dynamical system [15|. While it is very 
difficult to see what happens in the general situation, we know that the mean value of the 
Melnikov function is zero if the mean value of all the external perturbations is zero, as is 
shown by Eq. (153]) . provided it is well defined in this case. This suggests to us that in a 
broad number of situations the topology of the phase space will remain unchanged if the 
mean value of all the external perturbations vanishes, even when some of its fixed points 
were nonhyperbolic. It is, however, very difficult to prove a result in that direction, since 
perturbating a system in a critical state could have unpredictable consequences. Further- 
more, the existence of the Fourier series expansion of the Melnikov function is no longer 
guaranteed, as the loss of hyperbolicity implies that we do not have an exponential decay 
at infinity of the absolute value in Eq. fl56l) . 

If some of the perturbations have a nonzero mean, then the problem becomes much 
more difficult. The system will be moved out of criticality at some of its points, and its 
behavior will be more similar to that of the (partially) noncritical phase corresponding to 
the autonomous system with the parameters retuned according to the mean value of the 
external perturbations. Unfortunately, we cannot map this situation to that described in 
Sec. Ull because in this case we have no control on the amplitude of the perturbation, which 
may be comparable to the distance to some of the critical points. So the problem becomes 
genuinely nonperturbational, and we can no longer employ a technique like the Melnikov 
function. In this case it is difficult to deal even with very slow perturbations, since this type 
of settings favors the appearance of soft modes, which rules out the possibility of treating 
the external signals adiabatically 221 ]. 

Another problem we would like to deal with is the case of several perturbations with 
different periods. The simplest case is that in which the ratio of the periods of all possible 
different pairs of perturbations that are affecting the system is a rational number. In this 
case, we say that the perturbations are commensurable, and we can find a period which is 
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common to all of them. Once obtained, we have reduced our problem to the one studied 
in Sec [TTJ The case of incommensurable perturbations is, of course, more complex, as it 
implies that the dynamics of the perturbed system cannot be reduced to the dynamics of a 
periodic Poincare map. The fixed points of the unperturbed system give rise to quasiperiodic 
solutions of the perturbed problem, with their invariant manifolds being quasiperiodic as 

weii y . 

There are many other questions that remain to be answered. One is the possible appear- 
ance of chaotic behavior induced by internal stochasticity. As noted in [3], for two reacting 
species we have a four-dimensional eikonal Hamiltonian, which will give rise, in general, to 
a three-dimensional dynamical system on some Riemannian manifold. In this case, chaos, 
unlike in the two-dimensional mean-field dynamical system, is indeed possible. But chaos is 
also possible in the case of one reacting species obeying reaction rules with time- dependent 
rates. In this situation, "energy" is no longer conserved and the nonautonomous dynamical 
system, without integrals of motion, becomes effectively three dimensional. Indeed, situa- 
tions like the one plotted in Fig. (TSJ), that is, when the heteroclinic connection is preserved 
due to the intersection of the stable and unstable manifolds, give rise to complex dynami- 
cal scenarios. These include the appearance of Smale horseshoes and even the presence of 



strange attractors related to the creation and destruction of such horseshoes [30l . l3ll . [32]. 
Although the generation of chaotic behavior of this type can be attractive from a nonlinear 
dynamics point of view, we are not aware, at this point, of its possible physical meaning. 

Of course, studying multispecies reactions, both autonomous and nonautonomous, is a 
very interesting problem that deserves further efforts. The nonautonomous situation implies 
the technical difficulty of extending the Melnikov function theory to a higher dimensionality. 



33, 



34J ] could perhaps be adapted for 



The ideas developed for three-dimensional systems 
studying four- or higher- dimensional problems and to try to obtain in this way the results 
that we would need to understand the more general multispecies reactions. Also, as we have 
already pointed out, the Melnikov function is a perturbative result. It allows us to treat 
arbitrary frequency but necessarily small perturbations. To fully understand the dynamics 
of nonautonomous chemical reactions, even in the one species case, we would need a nonper- 
turbational result. With it we could try to address questions such as the forcing of systems 
near criticality or the appearance of soft modes. As we can see, chemical kinetics presents 
a very complex phenomenology, and it is challenging from a mathematical point of view. A 
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deep understanding of the physics of these nonequilibrium systems will presumably imply 
the parallel development of powerful methodological techniques. 
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APPENDIX: SECOND-ORDER MELNIKOV FUNCTION 

Consider the two-dimensional dynamical system 

x = f(x) + eg(x,t) + e 2 h(x,t), (A.l) 

where x = (xi,x 2 ), f = (A,/^), g = (91,92), and h = (hi,h 2 ) are sufficiently regular 
functions. The second-order Melnikov function reads 

/oo 
f(x fc (t)) A [2- 1 J D 2 f(x,(t)){xi(t + t ,t )} 2 
-00 

+£>g(x ft (t), t + £oK(t + t , t ) + h(x h (t), t + t )] dt, (A.2) 

where D 2 i and Dg denote the Hessian and Jacobian of f and g, respectively, x^(t) is the 
solution parametrizing some given heteroclinic connection, and x^ is the solution to the 
variational equation 

±l(t, t ) = Df(x h (t - t ))x£(i, to) + gfaft - t ), t). (A.3) 

We now assume that this dynamical system is a reaction Hamiltonian system and both 
g and h come from small nonautonomous perturbations of the Hamiltonian parameters, 
just like in Sec. HI1 Then, because Eq. (I A. 3 1) is linear, its solution will depend linearly on 
the Fourier modes of the time-dependent parameters evaluated at t . But this solution 
appears quadratically in Eq. flA.2j) and also multiplying the Jacobian of g, which depends 
linearly, by assumption, on the same Fourier modes. This quadratic dependence complicates 
recasting the second-order Melnikov function into a form similar to Eq. ( 15^1) . which would 
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help enormously for calculating its mean value. 
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